############
## Price of Peace Replication 1.R
## Last modified: April 22, 2019
## Josh Kertzer

## This file generates the analysis for the main text - run Price of Replication 0.R first to generate the necessary dataframes, etc.
## All of the following analyses were carried out using R version 3.5.2 GUI 1.7.0 El Capitan build (7612) on a 2017 iMac Pro running macOS Mojave 10.14.4

############


#First, just look at cell means

set.seed(43215)
B <- 1500

#### Figure 1 in paper: cell means

combined.trust.mat <- matrix(NA, nrow=B, ncol=9)
combined.danger.mat <- matrix(NA, nrow=B, ncol=9)
combined.support.mat <- matrix(NA, nrow=B, ncol=9)
combined.likely.mat <- matrix(NA, nrow=B, ncol=9)

for (j in 1:B){
	k <- sample(1:nrow(dat1),nrow(dat1),replace=TRUE)
	temp <- dat1[k,]
	combined.trust.mat[j,1] <- mean(temp$trustworthy[which(temp$Inspection==1 & temp$FuelCycle==1)], na.rm=TRUE)
	combined.trust.mat[j,2] <- mean(temp$trustworthy[which(temp$Inspection==2 & temp$FuelCycle==1)], na.rm=TRUE)
	combined.trust.mat[j,3] <- mean(temp$trustworthy[which(temp$Inspection==3 & temp$FuelCycle==1)], na.rm=TRUE)
	combined.trust.mat[j,4] <- mean(temp$trustworthy[which(temp$Inspection==1 & temp$FuelCycle==2)], na.rm=TRUE)
	combined.trust.mat[j,5] <- mean(temp$trustworthy[which(temp$Inspection==2 & temp$FuelCycle==2)], na.rm=TRUE)
	combined.trust.mat[j,6] <- mean(temp$trustworthy[which(temp$Inspection==3 & temp$FuelCycle==2)], na.rm=TRUE)	
	combined.trust.mat[j,7] <- mean(temp$trustworthy[which(temp$Inspection==1 & temp$FuelCycle==3)], na.rm=TRUE)
	combined.trust.mat[j,8] <- mean(temp$trustworthy[which(temp$Inspection==2 & temp$FuelCycle==3)], na.rm=TRUE)
	combined.trust.mat[j,9] <- mean(temp$trustworthy[which(temp$Inspection==3 & temp$FuelCycle==3)], na.rm=TRUE)
	combined.danger.mat[j,1] <- mean(temp$dangerous[which(temp$Inspection==1 & temp$FuelCycle==1)], na.rm=TRUE)
	combined.danger.mat[j,2] <- mean(temp$dangerous[which(temp$Inspection==2 & temp$FuelCycle==1)], na.rm=TRUE)
	combined.danger.mat[j,3] <- mean(temp$dangerous[which(temp$Inspection==3 & temp$FuelCycle==1)], na.rm=TRUE)
	combined.danger.mat[j,4] <- mean(temp$dangerous[which(temp$Inspection==1 & temp$FuelCycle==2)], na.rm=TRUE)
	combined.danger.mat[j,5] <- mean(temp$dangerous[which(temp$Inspection==2 & temp$FuelCycle==2)], na.rm=TRUE)
	combined.danger.mat[j,6] <- mean(temp$dangerous[which(temp$Inspection==3 & temp$FuelCycle==2)], na.rm=TRUE)	
	combined.danger.mat[j,7] <- mean(temp$dangerous[which(temp$Inspection==1 & temp$FuelCycle==3)], na.rm=TRUE)
	combined.danger.mat[j,8] <- mean(temp$dangerous[which(temp$Inspection==2 & temp$FuelCycle==3)], na.rm=TRUE)
	combined.danger.mat[j,9] <- mean(temp$dangerous[which(temp$Inspection==3 & temp$FuelCycle==3)], na.rm=TRUE)
	combined.support.mat[j,1] <- mean(temp$dealSupport[which(temp$Inspection==1 & temp$FuelCycle==1)], na.rm=TRUE)
	combined.support.mat[j,2] <- mean(temp$dealSupport[which(temp$Inspection==2 & temp$FuelCycle==1)], na.rm=TRUE)
	combined.support.mat[j,3] <- mean(temp$dealSupport[which(temp$Inspection==3 & temp$FuelCycle==1)], na.rm=TRUE)
	combined.support.mat[j,4] <- mean(temp$dealSupport[which(temp$Inspection==1 & temp$FuelCycle==2)], na.rm=TRUE)
	combined.support.mat[j,5] <- mean(temp$dealSupport[which(temp$Inspection==2 & temp$FuelCycle==2)], na.rm=TRUE)
	combined.support.mat[j,6] <- mean(temp$dealSupport[which(temp$Inspection==3 & temp$FuelCycle==2)], na.rm=TRUE)	
	combined.support.mat[j,7] <- mean(temp$dealSupport[which(temp$Inspection==1 & temp$FuelCycle==3)], na.rm=TRUE)
	combined.support.mat[j,8] <- mean(temp$dealSupport[which(temp$Inspection==2 & temp$FuelCycle==3)], na.rm=TRUE)
	combined.support.mat[j,9] <- mean(temp$dealSupport[which(temp$Inspection==3 & temp$FuelCycle==3)], na.rm=TRUE)	
	combined.likely.mat[j,1] <- mean(temp$dealLikely[which(temp$Inspection==1 & temp$FuelCycle==1)], na.rm=TRUE)
	combined.likely.mat[j,2] <- mean(temp$dealLikely[which(temp$Inspection==2 & temp$FuelCycle==1)], na.rm=TRUE)
	combined.likely.mat[j,3] <- mean(temp$dealLikely[which(temp$Inspection==3 & temp$FuelCycle==1)], na.rm=TRUE)
	combined.likely.mat[j,4] <- mean(temp$dealLikely[which(temp$Inspection==1 & temp$FuelCycle==2)], na.rm=TRUE)
	combined.likely.mat[j,5] <- mean(temp$dealLikely[which(temp$Inspection==2 & temp$FuelCycle==2)], na.rm=TRUE)
	combined.likely.mat[j,6] <- mean(temp$dealLikely[which(temp$Inspection==3 & temp$FuelCycle==2)], na.rm=TRUE)	
	combined.likely.mat[j,7] <- mean(temp$dealLikely[which(temp$Inspection==1 & temp$FuelCycle==3)], na.rm=TRUE)
	combined.likely.mat[j,8] <- mean(temp$dealLikely[which(temp$Inspection==2 & temp$FuelCycle==3)], na.rm=TRUE)
	combined.likely.mat[j,9] <- mean(temp$dealLikely[which(temp$Inspection==3 & temp$FuelCycle==3)], na.rm=TRUE)
	}

combined.mat <- data.frame(cbind(rbind(cbind(apply(combined.trust.mat, 2,mean),apply(combined.trust.mat, 2,quantile,0.025), apply(combined.trust.mat, 2,quantile,0.05), apply(combined.trust.mat, 2,quantile,0.95), apply(combined.trust.mat, 2,quantile,0.975)), cbind(apply(combined.danger.mat, 2,mean),apply(combined.danger.mat, 2,quantile,0.025), apply(combined.danger.mat, 2,quantile,0.05), apply(combined.danger.mat, 2,quantile,0.95), apply(combined.danger.mat, 2,quantile,0.975)), cbind(apply(combined.support.mat, 2,mean),apply(combined.support.mat, 2,quantile,0.025), apply(combined.support.mat, 2,quantile,0.05), apply(combined.support.mat, 2,quantile,0.95), apply(combined.support.mat, 2,quantile,0.975)), cbind(apply(combined.likely.mat, 2,mean),apply(combined.likely.mat, 2,quantile,0.025), apply(combined.likely.mat, 2,quantile,0.05), apply(combined.likely.mat, 2,quantile,0.95), apply(combined.likely.mat, 2,quantile,0.975))), Inspections=rep(c("1. None", "2. Civilian Only", "3. All Sites"),4), Fuel.Cycle=rep(c("1. Continue", "2. Freeze", "3. Dismantle"),each=3), DV=rep(c("A. Trustworthiness of Iran", "B. Danger posed by Iran", "C. Support for deal", "D. Likelihood of Iran making deal"), each=9)))
colnames(combined.mat) <- c("Y", "Low_95", "Low_90", "Hi_90", "Hi_95", "Inspections","Fuel.Cycle", "DV")  
combined.mat$Y <- as.numeric(as.character(combined.mat$Y))
combined.mat$Low_95 <- as.numeric(as.character(combined.mat$Low_95))
combined.mat$Hi_95 <- as.numeric(as.character(combined.mat$Hi_95))
combined.mat$Low_90 <- as.numeric(as.character(combined.mat$Low_90))
combined.mat$Hi_90 <- as.numeric(as.character(combined.mat$Hi_90))

dev.new(height=7,width=8)
qplot(Inspections, Y, data=combined.mat, color=Fuel.Cycle, ymin=Low_95, ymax=Hi_95, ylab="", xlab="Inspection Regime", geom="pointrange") + facet_wrap(~DV) + theme_bw() + geom_line(aes(x=Inspections,y=Y, color=Fuel.Cycle, group=Fuel.Cycle)) + scale_color_discrete(name="Fuel Cycle")


##### Table 1: Trustworthiness of Iran

mod.t1 <- lm(trustworthy ~ Inspection_NotAllow + Inspection_All + FuelCycle_Continue + FuelCycle_Dismantle, data=dat1)
mod.t2 <- lm(trustworthy ~ (Inspection_NotAllow + Inspection_All) + (FuelCycle_Continue + FuelCycle_Dismantle) + mi1 + ci1 + iso1 + Male + log(age) + White + educ + ideo1, data=dat1)
mod.t3 <- lm(trustworthy ~ (Inspection_NotAllow + Inspection_All) + (FuelCycle_Continue + FuelCycle_Dismantle) + mi1 + ci1 + iso1 + iran + Male + log(age) + White + educ + ideo1, data=dat1)
mod.t4 <- lm(trustworthy ~ (Inspection_NotAllow + Inspection_All) + (FuelCycle_Continue + FuelCycle_Dismantle) + mi1 + ci1 + iso1 + iran + cognition1 + Male + log(age) + White + educ + ideo1, data=dat1)

stargazer(mod.t1, mod.t2, mod.t3, mod.t4,  title="Trustworthiness of Iran", omit.stat=c("LL", "ser", "f"), style="apsr", digits=3, label="tab:trust1", covariate.labels=c("Inspections: None", "Inspections: All", "Fuel Cycle: Continue", "Fuel Cycle: Dismantle", "Militant internationalism", "Cooperative internationalism", "Isolationism", "Feelings towards Iran", "Need for cognition", "Male", "log(Age)", "White", "Education", "Ideology", "Constant"))

##### Table 2: CI moderates the impact of costly signals

mod.ci1 <- lm(trustworthy ~ (Inspection_NotAllow + Inspection_All)*ci1 + (FuelCycle_Continue + FuelCycle_Dismantle)*ci1 + mi1 + iso1 + iran + cognition1 + Male + log(age) + White + educ + ideo1, data=dat1)
mod.ci2 <- lm(dangerous ~ (Inspection_NotAllow + Inspection_All)*ci1 + (FuelCycle_Continue + FuelCycle_Dismantle)*ci1 + mi1 + iso1 + iran + cognition1+ Male + log(age) + White + educ + ideo1, data=dat1)
mod.ci3 <- lm(dealSupport ~ (Inspection_NotAllow + Inspection_All)*ci1 + (FuelCycle_Continue + FuelCycle_Dismantle)*ci1 + mi1 + iso1 + iran + cognition1 + Male + log(age) + White + educ + ideo1, data=dat1)

stargazer(mod.ci1, mod.ci2, mod.ci3,  title="CI moderates the impact of costly signals", omit.stat=c("LL", "ser", "f"), style="apsr", digits=3, label="tab:ci1", covariate.labels=c("Inspections: None", "Inspections: All", "Cooperative internationalism", "Fuel Cycle: Continue", "Fuel Cycle: Dismantle", "Militant internationalism", "Isolationism", "Feelings towards Iran", "Need for cognition", "Male", "log(Age)", "White", "Education", "Ideology", "Inspections: None x CI", "Inspections: All x CI", "Fuel Cycle: Continue x CI", "Fuel Cycle: Dismantle x CI", "Constant"))

##### Figure 3: Conditional effects of signals by cooperative internationalism

set.seed(43215)
B <- 1500

ci.mat1a <- matrix(NA, nrow=B,ncol=2)
for (j in 1:B){
	k <- sample(1:nrow(dat1),nrow(dat1),replace=TRUE)
	temp <- dat1[k,]
	temp.mod <- lm(trustworthy ~ (Inspection_NoMilitary + Inspection_All)*ci1 + (FuelCycle_Freeze + FuelCycle_Dismantle) + mi1 + iso1 + Male + log(age) + White + educ + ideo1, data=temp)
	ci.mat1a[j,] <- coef(temp.mod)[c(3,15)]
}

ci.mat1b <- matrix(NA, nrow=B,ncol=2)
for (j in 1:B){
	k <- sample(1:nrow(dat1),nrow(dat1),replace=TRUE)
	temp <- dat1[k,]
	temp.mod <- lm(trustworthy ~ (Inspection_NoMilitary + Inspection_All) + (FuelCycle_Freeze + FuelCycle_Dismantle)*ci1 + mi1 + iso1 + Male + log(age) + White + educ + ideo1, data=temp)
	ci.mat1b[j,] <- coef(temp.mod)[c(5,15)]
}

ci.mat2a <- matrix(NA, nrow=B,ncol=2)
for (j in 1:B){
	k <- sample(1:nrow(dat1),nrow(dat1),replace=TRUE)
	temp <- dat1[k,]
	temp.mod <- lm(dangerous ~ (Inspection_NoMilitary + Inspection_All)*ci1 + (FuelCycle_Freeze + FuelCycle_Dismantle) + mi1 + iso1 + Male + log(age) + White + educ + ideo1, data=temp)
	ci.mat2a[j,] <- coef(temp.mod)[c(3,15)]
}

ci.mat2b <- matrix(NA, nrow=B,ncol=2)
for (j in 1:B){
	k <- sample(1:nrow(dat1),nrow(dat1),replace=TRUE)
	temp <- dat1[k,]
	temp.mod <- lm(dangerous ~ (Inspection_NoMilitary + Inspection_All) + (FuelCycle_Freeze + FuelCycle_Dismantle)*ci1 + mi1 + iso1 + Male + log(age) + White + educ + ideo1, data=temp)
	ci.mat2b[j,] <- coef(temp.mod)[c(5,15)]
}

ci.mat3a <- matrix(NA, nrow=B,ncol=2)
for (j in 1:B){
	k <- sample(1:nrow(dat1),nrow(dat1),replace=TRUE)
	temp <- dat1[k,]
	temp.mod <- lm(dealSupport ~ (Inspection_NoMilitary + Inspection_All)*ci1 + (FuelCycle_Freeze + FuelCycle_Dismantle) + mi1 + iso1 + Male + log(age) + White + educ + ideo1, data=temp)
	ci.mat3a[j,] <- coef(temp.mod)[c(3,15)]
}

ci.mat3b <- matrix(NA, nrow=B,ncol=2)
for (j in 1:B){
	k <- sample(1:nrow(dat1),nrow(dat1),replace=TRUE)
	temp <- dat1[k,]
	temp.mod <- lm(dealSupport ~ (Inspection_NoMilitary + Inspection_All) + (FuelCycle_Freeze + FuelCycle_Dismantle)*ci1 + mi1 + iso1 + Male + log(age) + White + educ + ideo1, data=temp)
	ci.mat3b[j,] <- coef(temp.mod)[c(5,15)]
}

c.val <- seq(0,1,length.out=11)
ci.plot1a <- matrix(NA, 11,3)
for (i in 1:11){
	ci.plot1a[i,1] <- mean(ci.mat1a[,1] + c.val[i]*ci.mat1a[,2])
	ci.plot1a[i,2] <- quantile(ci.mat1a[,1] + c.val[i]*ci.mat1a[,2], 0.025)
	ci.plot1a[i,3] <- quantile(ci.mat1a[,1] + c.val[i]*ci.mat1a[,2], 0.975)  
}

ci.plot1b <- matrix(NA, 11,3)
for (i in 1:11){
	ci.plot1b[i,1] <- mean(ci.mat1b[,1] + c.val[i]*ci.mat1b[,2])
	ci.plot1b[i,2] <- quantile(ci.mat1b[,1] + c.val[i]*ci.mat1b[,2], 0.025)
	ci.plot1b[i,3] <- quantile(ci.mat1b[,1] + c.val[i]*ci.mat1b[,2], 0.975)  
}

ci.plot2a <- matrix(NA, 11,3)
for (i in 1:11){
	ci.plot2a[i,1] <- mean(ci.mat2a[,1] + c.val[i]*ci.mat2a[,2])
	ci.plot2a[i,2] <- quantile(ci.mat2a[,1] + c.val[i]*ci.mat2a[,2], 0.025)
	ci.plot2a[i,3] <- quantile(ci.mat2a[,1] + c.val[i]*ci.mat2a[,2], 0.975)  
}

ci.plot2b <- matrix(NA, 11,3)
for (i in 1:11){
	ci.plot2b[i,1] <- mean(ci.mat2b[,1] + c.val[i]*ci.mat2b[,2])
	ci.plot2b[i,2] <- quantile(ci.mat2b[,1] + c.val[i]*ci.mat2b[,2], 0.025)
	ci.plot2b[i,3] <- quantile(ci.mat2b[,1] + c.val[i]*ci.mat2b[,2], 0.975)  
}

ci.plot3a <- matrix(NA, 11,3)
for (i in 1:11){
	ci.plot3a[i,1] <- mean(ci.mat3a[,1] + c.val[i]*ci.mat3a[,2])
	ci.plot3a[i,2] <- quantile(ci.mat3a[,1] + c.val[i]*ci.mat3a[,2], 0.025)
	ci.plot3a[i,3] <- quantile(ci.mat3a[,1] + c.val[i]*ci.mat3a[,2], 0.975)  
}

ci.plot3b <- matrix(NA, 11,3)
for (i in 1:11){
	ci.plot3b[i,1] <- mean(ci.mat3b[,1] + c.val[i]*ci.mat3b[,2])
	ci.plot3b[i,2] <- quantile(ci.mat3b[,1] + c.val[i]*ci.mat3b[,2], 0.025)
	ci.plot3b[i,3] <- quantile(ci.mat3b[,1] + c.val[i]*ci.mat3b[,2], 0.975)  
}

#Create separate plots for each panel
dev.new(height=4,width=4)
par(mar=c(4.1,4.1,1.1,1.1), oma=c(0,0,0,0))

plot(c.val, ci.plot1a[,1], type="n", xlab="Cooperative internationalism", ylab="Conditional effect of complete inspections", main="", ylim=c(-2,4))
polygon(c(c.val, rev(c.val)), c(ci.plot1a[,2], rev(ci.plot1a[,3])), col="grey")
lines(c.val, ci.plot1a[,1], lwd=2)
abline(h=0, lty=3)

plot(c.val, ci.plot1b[,1], type="n", xlab="Cooperative internationalism", ylab="Conditional effect of dismantling fuel cycle", main=" ", ylim=c(-2,4))
polygon(c(c.val, rev(c.val)), c(ci.plot1b[,2], rev(ci.plot1b[,3])), col="grey")
lines(c.val, ci.plot1b[,1], lwd=2)
abline(h=0, lty=3)

plot(c.val, ci.plot2a[,1], type="n", xlab="Cooperative internationalism", ylab="Conditional effect of complete inspections", main="", ylim=c(-3,2))
polygon(c(c.val, rev(c.val)), c(ci.plot2a[,2], rev(ci.plot2a[,3])), col="grey")
lines(c.val, ci.plot2a[,1], lwd=2)
abline(h=0, lty=3)

plot(c.val, ci.plot2b[,1], type="n", xlab="Cooperative internationalism", ylab="Conditional effect of dismantling fuel cycle", main="", ylim=c(-3,2))
polygon(c(c.val, rev(c.val)), c(ci.plot2b[,2], rev(ci.plot2b[,3])), col="grey")
lines(c.val, ci.plot2b[,1], lwd=2)
abline(h=0, lty=3)

plot(c.val, ci.plot3a[,1], type="n", xlab="Cooperative internationalism", ylab="Conditional effect of complete inspections", main="", ylim=c(-1,3))
polygon(c(c.val, rev(c.val)), c(ci.plot3a[,2], rev(ci.plot3a[,3])), col="grey")
lines(c.val, ci.plot3a[,1], lwd=2)
abline(h=0, lty=3)

plot(c.val, ci.plot3b[,1], type="n", xlab="Cooperative internationalism", ylab="Conditional effect of dismantling fuel cycle", main="", ylim=c(-1,3))
polygon(c(c.val, rev(c.val)), c(ci.plot3b[,2], rev(ci.plot3b[,3])), col="grey")
lines(c.val, ci.plot3b[,1], lwd=2)
abline(h=0, lty=3)

##### Figure 4: Within-subject effects of signals on Iranian trustworthiness

# Calculate fitted values from each model
set.seed(43215)
B <- 1500

dat2a <- subset(dat2, dat2$complied==1)
dat2b <- subset(dat2, dat2$ballistic==1)
dat2c <- subset(dat2, dat2$pledge==1)

x <- seq(0,1,length.out=6) #Values of CI
v <- c(1,0,mean(dat2$mi1, na.rm=TRUE), mean(dat2$iran1, na.rm=TRUE), 0.5) #Values to hold other covariates at

boot.out1 <- matrix(NA, nrow=B, ncol=6)
for (i in 1:nrow(boot.out1)){
	j <- sample(1:nrow(dat2), nrow(dat2), replace=TRUE)
	temp <- dat2[j,]
	mod.1a1 <- lm(trustDiff ~ ci1 + mi1 + iran1 + news_iran1, data=temp)
	temp.b <- coef(mod.1a1)
	boot.out1[i,] <- as.vector(temp.b %*% v) + (temp.b[2]*x) #Predicted value
	}

boot.out2 <- matrix(NA, nrow=B, ncol=6)
for (i in 1:nrow(boot.out2)){
	j <- sample(1:nrow(dat2a), nrow(dat2a), replace=TRUE)
	temp <- dat2a[j,]
	mod.2a1 <- lm(trustDiff ~ ci1 + mi1 + iran1 + news_iran1, data=temp)
	temp.b <- coef(mod.2a1)
	boot.out2[i,] <- as.vector(temp.b %*% v) + (temp.b[2]*x) #Predicted value
	}
	
boot.out3 <- matrix(NA, nrow=B, ncol=6)
for (i in 1:nrow(boot.out3)){
	j <- sample(1:nrow(dat2b), nrow(dat2b), replace=TRUE)
	temp <- dat2b[j,]
	mod.2b1 <- lm(trustDiff ~ ci1 + mi1 + iran1 + news_iran1, data=temp)
	temp.b <- coef(mod.2b1)
	boot.out3[i,] <- as.vector(temp.b %*% v) + (temp.b[2]*x) #Predicted value
	}
	
boot.out4 <- matrix(NA, nrow=B, ncol=6)
for (i in 1:nrow(boot.out4)){
	j <- sample(1:nrow(dat2c), nrow(dat2c), replace=TRUE)
	temp <- dat2c[j,]
	mod.2c1 <- lm(trustDiff ~ ci1 + mi1 + iran1 + news_iran1, data=temp)
	temp.b <- coef(mod.2c1)
	boot.out4[i,] <- as.vector(temp.b %*% v) + (temp.b[2]*x) #Predicted value
	}
		
out.mat <- cbind(apply(boot.out1,2,mean), apply(boot.out1,2,quantile,0.025), apply(boot.out1,2,quantile,0.975))
out.mat2 <- cbind(apply(boot.out2,2,mean), apply(boot.out2,2,quantile,0.025), apply(boot.out2,2,quantile,0.975))
out.mat3 <- cbind(apply(boot.out3,2,mean), apply(boot.out3,2,quantile,0.025), apply(boot.out3,2,quantile,0.975))
out.mat4 <- cbind(apply(boot.out4,2,mean), apply(boot.out4,2,quantile,0.025), apply(boot.out4,2,quantile,0.975))

combined.mat <- as.data.frame(rbind(out.mat, out.mat2, out.mat3, out.mat4))
colnames(combined.mat) <- c("y", "low", "high")
combined.mat$x <- rep(x,4)
combined.mat$DV <- rep(c("1. Pooled Model", "2. Compliance Condition", "3. Ballistic Condition", "4. Pledge Condition"),each=6)

dev.new(height=5,width=5)
ggplot(aes(x, y), data=combined.mat) + geom_line() + facet_wrap(~DV) + geom_ribbon(aes(ymin=low, ymax=high), alpha=I(0.3)) + theme_bw() + labs(x="Cooperative Internationalism", y="Within-subject updating of Iranian trustworthiness") + geom_hline(yintercept=0, lty=3)

###### Figure 5: Costly signals induce polarization

#Panel a: 2015 experiment

B <- 1500
set.seed(43215)

var.arr <- array(NA, dim=c(B,3,3))
for (i in 1:B){
	j <- sample(seq_len(nrow(dat1)),nrow(dat1), replace=TRUE)
	temp <- dat1[j,]
	var.arr[i,1,1] <- var(temp$trustworthy[which(temp$FuelCycle_Continue==1 & temp$Inspection_NotAllow==1)], na.rm=TRUE)
	var.arr[i,2,1] <- var(temp$trustworthy[which(temp$FuelCycle_Freeze==1 & temp$Inspection_NotAllow==1)], na.rm=TRUE)
	var.arr[i,3,1] <- var(temp$trustworthy[which(temp$FuelCycle_Dismantle==1 & temp$Inspection_NotAllow==1)], na.rm=TRUE)
	var.arr[i,1,2] <- var(temp$trustworthy[which(temp$FuelCycle_Continue==1 & temp$Inspection_NoMilitary==1)], na.rm=TRUE)
	var.arr[i,2,2] <- var(temp$trustworthy[which(temp$FuelCycle_Freeze==1 & temp$Inspection_NoMilitary==1)], na.rm=TRUE)
	var.arr[i,3,2] <- var(temp$trustworthy[which(temp$FuelCycle_Dismantle==1 & temp$Inspection_NoMilitary==1)], na.rm=TRUE)
	var.arr[i,1,3] <- var(temp$trustworthy[which(temp$FuelCycle_Continue==1 & temp$Inspection_All==1)], na.rm=TRUE)
	var.arr[i,2,3] <- var(temp$trustworthy[which(temp$FuelCycle_Freeze==1 & temp$Inspection_All==1)], na.rm=TRUE)
	var.arr[i,3,3] <- var(temp$trustworthy[which(temp$FuelCycle_Dismantle==1 & temp$Inspection_All==1)], na.rm=TRUE)	
}

dimnames(var.arr)[[2]] <- c("Fuel Cycle: Continue", "Fuel Cycle: Freeze", "Fuel Cycle: Dismantle")
dimnames(var.arr)[[3]] <- c("1. None", "2. Civilian Only", "3. All Sites")
var.arr2 <- melt(var.arr)
colnames(var.arr2) <- c("Var1", "Fuel.Cycle", "Inspections", "value")

dev.new(height=3.5,width=8)
ggplot(var.arr2, aes(y=value, x=Inspections)) + geom_violin(fill="grey") + facet_wrap(~Fuel.Cycle, ncol=3) + theme_bw() + labs(y="Variance estimates of trustworthiness", x="Inspection Regime")

#Panel b: 2018 experiment

B <- 1500
set.seed(43215)

var.mat <- matrix(NA, B, ncol=6)
for (i in 1:B){
	j <- sample(seq_len(nrow(dat2)),nrow(dat2), replace=TRUE)
	temp <- dat2[j,]
	var.mat[i,1] <- var(temp$trustworthyPrior[which(temp$complied==1)], na.rm=TRUE)
	var.mat[i,2] <- var(temp$trustworthy[which(temp$complied==1)], na.rm=TRUE)
	var.mat[i,3] <- var(temp$trustworthyPrior[which(temp$ballistic==1)], na.rm=TRUE)
	var.mat[i,4] <- var(temp$trustworthy[which(temp$ballistic==1)], na.rm=TRUE)
	var.mat[i,5] <- var(temp$trustworthyPrior[which(temp$pledge==1)], na.rm=TRUE)	
	var.mat[i,6] <- var(temp$trustworthy[which(temp$pledge==1)], na.rm=TRUE)				
}

colnames(var.mat) <- c("Compliance (Prior)", "Compliance (Posterior)", "Ballistic (Prior)", "Ballistic (Posterior)", "Pledge (Prior)", "Pledge (Posterior)")
var.mat2 <- melt(var.mat)
colnames(var.mat2) <- c("i","Signal", "value")
var.mat2$Prior <- as.factor(rep(c("Prior", "Posterior"), each=B))
var.mat2$Signal2 <- as.factor(rep(c("Compliance", "Ballistic", "Pledge"), each=2*B))
var.mat2$Prior <- relevel(var.mat2$Prior, ref="Prior")
var.mat2$Signal2 <- relevel(var.mat2$Signal2, ref="Compliance")

dev.new(height=3.5, width=8)
ggplot(var.mat2, aes(y=value, x=Prior)) + geom_violin(fill="grey") + facet_wrap(~Signal2) + theme_bw() + labs(y="Variance estimates of trustworthiness", x="")

## Bootstrapped test of difference 
#Study 1:
1 - length(which(var.arr[,1,1] < var.arr[,1,3]))/nrow(var.arr) #p < 0.047: All sites vs none | fuel cycle continue
1 - length(which(var.arr[,2,1] < var.arr[,2,3]))/nrow(var.arr) #p < 0.009: All sites vs none | fuel cycle freeze
1 - length(which(var.arr[,3,1] < var.arr[,3,3]))/nrow(var.arr) #p < 0.000: All sites vs none | fuel cycle dismantle

#Study 2: 
1 - length(which(var.mat[,1] < var.mat[,2]))/nrow(var.mat) #p < 0.000: None vs compliance
1 - length(which(var.mat[,3] < var.mat[,4]))/nrow(var.mat) #p < 0.015: None vs ballistic
1 - length(which(var.mat[,5] < var.mat[,6]))/nrow(var.mat) #p < 0.000: None vs pledge

